############## 
############## 
############## Models
remove(list = ls())

base::library(conflicted)
base::library(tidyverse)
conflict_prefer("filter","dplyr")
base::library(ggplot2)
base::library(dplyr)
base::library(here)
conflict_prefer("here", "here")
base::library(systemfit)
base::library(stargazer)
base::library(DataCombine)
base::library(haven)
base::library(nnet)
base::library(ggeffects)

data <- readRDS(here("Data", "cumulative_2006-2020.Rds"))
glimpse(data)
table(data$year)

data2 <- data %>% filter(year == 2006  | year == 2008 | year == 2010 |year == 2012 | 
                           year == 2014 |  year == 2016 | year == 2018 | year == 2020)

data2 <- data2 %>% 
  mutate(fips = as.numeric(county_fips)) %>% 
  filter(!is.na(fips))
summary(data2$fips)

load(here("Data", "Deaths_AllYears.RData"))
glimpse(drug_data)
load(here("Data", "Model_Data.RData"))
glimpse(final_data)


county_data <- final_data %>% select(fips, rural_urban, tot_pop, 
                                     unemp_rate, defl_pcincome) %>% 
  distinct(fips, .keep_all = T)

drug_data <- left_join(drug_data, county_data, by = "fips")

county_data <- drug_data %>% 
  select(year, fips, rural_urban, death_rate, unemp_rate, defl_pcincome, tot_pop,
         lag_death_rate)

ces_data <- left_join(data2, county_data, by = c("year", "fips"))
glimpse(ces_data)

table(ces_data$pid7, useNA = "always")

ces_data <- ces_data %>% 
  mutate(pid7 = ifelse(pid7 == 8 | is.na(pid7), 4, pid7)) %>% 
  mutate(pid7 = as.factor(pid7))
table(ces_data$voted_rep, useNA = "always")

ces_data <- ces_data %>% 
  mutate(outcome_choice = ifelse(voted_rep == "[Democrat / Candidate 1]", "Dem", "Abst"),
         outcome_choice = ifelse(voted_rep == "[Republican / Candidate 2]", "Rep", outcome_choice),
         outcome_choice = ifelse(is.na(voted_rep), NA, outcome_choice)) 
table(ces_data$outcome_choice, useNA = "always")

table(ces_data$year)
glimpse(ces_data)


ces_data <- ces_data %>% 
  mutate(party_id = ifelse(pid3 == 1, "Dem", NA),
         party_id = ifelse(pid3 == 2, "Rep", party_id),
         party_id = ifelse(pid3 == 3 | pid3 == 4 |
                             pid3 == 5, "Ind", party_id)) %>% 
  mutate(party_id = as.factor(party_id))


ces_data <- ces_data %>% 
  mutate(white = ifelse(race == 1, 1, 0)) %>% 
  mutate(black = ifelse(race == 2, 1, 0)) %>% 
  mutate(hispanic = ifelse(race == 3, 1, 0)) %>% 
  mutate(asian = ifelse(race == 4, 1, 0)) %>% 
  mutate(other_race = ifelse(race == 5 |
                               race == 6 |
                               race == 7 |
                               race == 8, 1, 0))

ces_data <- ces_data %>% 
  mutate(college = ifelse(educ == 3 |
                            educ == 4 | 
                            educ == 5 |
                            educ == 6, 1, 0))

ces_data <- ces_data %>% 
  mutate(income = ifelse(faminc == "Less than 10k", 1, NA),
         income = ifelse(faminc == "10k - 20k", 2, income),
         income = ifelse(faminc == "20k - 30k", 3, income),
         income = ifelse(faminc == "30k - 40k", 4, income),
         income = ifelse(faminc == "40k - 50k", 5, income),
         income = ifelse(faminc == "50k - 60k", 6, income),
         income = ifelse(faminc == "60k - 70k", 7, income),
         income = ifelse(faminc == "70k - 80k", 8, income),
         income = ifelse(faminc == "80k - 100k", 9, income),
         income = ifelse(faminc == "100k - 120k", 10, income),
         income = ifelse(faminc == "120k - 150k", 11, income),
         income = ifelse(faminc == "150k+", 12, income))

ces_data <- ces_data %>% 
  mutate(gender = as.factor(gender)) %>% 
  mutate(race = as.factor(race)) %>% 
  mutate(state = as.factor(state)) %>% 
  mutate(log_age = log(age)) %>% 
  mutate(year_elec = as.factor(year)) %>% 
  mutate(time_trend = ifelse(year == 2006, 1, NA),
         time_trend = ifelse(year == 2008, 2, time_trend),
         time_trend = ifelse(year == 2010, 3, time_trend),
         time_trend = ifelse(year == 2012, 4, time_trend),
         time_trend = ifelse(year == 2014, 5, time_trend),
         time_trend = ifelse(year == 2016, 6, time_trend),
         time_trend = ifelse(year == 2018, 7, time_trend),
         time_trend = ifelse(year == 2020, 8, time_trend)) %>% 
  mutate(midterm = ifelse(year == 2006 | year == 2010 | year == 2014 |
                           year == 2018, 1, 0))


ces_data <- ces_data %>% 
  mutate(log_death = log1p(death_rate)) %>% 
  mutate(log_pop = log(tot_pop))

rep_data <- ces_data %>% 
  filter(!is.na(outcome_choice)) %>% 
  mutate(rep_vote = ifelse(outcome_choice == "Rep", 1, 0))  %>% 
  mutate(rural_urban = as.factor(rural_urban))


table(ces_data$outcome_choice)
ces_data <- ces_data %>% 
  mutate(fac_outcome = factor(outcome_choice,
                              levels = c("Rep", "Dem", "Abst"),
                              labels = c("Rep", "Dem", "Abst"))) %>% 
  mutate(rural_urban = as.factor(rural_urban))


################################################################################
rep_data <- rep_data %>%
  mutate(party_id3 = ifelse(pid7 == 1 | pid7 == 2 | pid7 == 3, 
                            "Democrat", NA),
         party_id3 = ifelse(pid7 == 4, "Independent", party_id3),
         party_id3 = ifelse(pid7 == 5 | pid7 == 6 | pid7 == 7, 
                            "Republican", party_id3)) %>% 
  mutate(party_id3 = factor(party_id3,
                            levels = c("Democrat", "Independent",
                                       "Republican"))) %>% 
  mutate(election = as.factor(year))

glimpse(rep_data)
summary(rep_pid3 <- glm(rep_vote ~ death_rate*party_id3 + unemp_rate + defl_pcincome +
                          log_pop + rural_urban + defl_pcincome + 
                          gender + black + hispanic + asian + other_race + 
                          log_age + college + income + state + election,
                        family = binomial(link = "logit"), weights = weight_cumulative,
                        data = rep_data))


model_pid3 <- rep_pid3 %>% 
  ggeffect(terms = c("death_rate[0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120]", 
                     "party_id3"), 
           ci.lvl = 0.90)

glimpse(model_pid3)
#table(model_drugs$response.level)
class(model_pid3)

model_pid3 <- data.frame(model_pid3) 


plot_partyid3 <- model_pid3 %>%
  mutate(group = factor(group, levels = c("Democrat",
                                          "Independent",
                                          "Republican"))) %>% 
  ggplot(aes(color = group, shape = group)) + 
  geom_pointrange(aes(x = x, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85) + 
  scale_shape_manual(name = "Party Identification:",
                     labels = c("Democrat",
                                "Independent",
                                "Republican"),
                     values = c(15, 16, 17)) +
  scale_color_manual(name = "Party Identification:",
                     labels = c("Democrat",
                                "Independent",
                                "Republican"),
                     values = c("blue",  "gray40", "red")) +
  theme_light() + 
  labs(x = "Overdose Death Rate", y = "Pred. Prob. of Republican Vote",
       title = NULL) +
  scale_x_continuous(breaks = c(0, 20, 40, 60, 80, 100, 120), lim = c(-1, 121)) +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), lim = c(0, 1.01)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 11),
        axis.text.y = element_text(size = 11),
        axis.title.y = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = "bottom",
        legend.direction = "horizontal",
        plot.title = element_text(size = 15))
plot_partyid3

summary(rep_pid3)
stargazer(rep_pid3,
          type = "latex",
          title = "",
          dep.var.labels = c("Prob. Voting for a Republican House Candidate"),
          no.space = FALSE, single.row = T,
          covariate.labels = c("Overdose Death Rate",
                               "PID3 - Independent",
                               "PID3 - Republican",
                               "Local Unemployment Rate",
                               "Per Capital Income (County",
                               "Total Population (Log)",
                               "Gender - Woman",
                               "Race - Black",
                               "Race - Hispanic",
                               "Race - Asian",
                               "Race - Other",
                               "Age (Log)",
                               "College Educated",
                               "Household Income",
                               "Death Rate X PID3 - Independent",
                               "Death Rate X PID3 - Republican",
                               "Constant"),
          omit = c("rural_urban", "state", "election"),
          omit.labels = c("Rural-Urban Code",
                          "State Fixed Effects",
                          "Election Fixed Effects"))

